#ifndef ARITHMETIC_H
#define ARITHMETIC_H

#include <boost/numeric/interval.hpp>
#include <fstream>
#include <vector>
#include <utility>


using namespace boost::numeric;
typedef interval<double> interval_d;
double eps_m = std::numeric_limits<double>::epsilon();

class RandomInput
{
   int seed;
   static const int a = 1664525;
   static const int c = 1013904223;
public:
   RandomInput(int seed) :seed(seed)
   { }

   // returns random integer in range 0 .. 2^31-1
   int get_int()
   {
      return (seed = seed * a + c) & 0x7fffffff;
   }

   double get_double()
   {
      const double lo = -1000;
      const double hi = 1000;
      const double scale = 1. / (1 << 30);
      return (lo + (hi - lo) * scale * (get_int() & 0x3fffffff));
   }
};


struct point
{
   double x;
   double y;
   point() { };
   point (double x, double y): x(x), y(y) { }
   point(RandomInput& rand):  x(rand.get_double()), y(rand.get_double()) { }
};

std::istream& operator>>(std::istream& in, point& a)
{
   in >> a.x >> a.y;
   return in;
}

//adaptive arithmetic
std::pair<double, double> sum(double a, double b)
{
   double res = a + b, roundoff;
   double bv = res - a, av = res - bv;
   double br = b - bv, ar = a - av;
   roundoff = ar + br;

   return std::make_pair(res, roundoff);
}

std::pair<double, double> split(double a, int s)
{
   double c = ((1LL << s) + 1LL) * a, ab = c - a;
   double hi = c - ab, lo = a - hi;

   return std::make_pair(hi, lo);
}

std::pair<double, double> multiply(double a, double b)
{
   double res = a * b, 
      roundoff;

   int s = std::numeric_limits<double>:: digits / 2
      + std::numeric_limits<double>:: digits % 2 ;

   double a_hi, a_lo, b_hi, b_lo;

   a_hi = split(a, s).first, a_lo = split(a, s).second;
   b_hi = split(b, s).first, b_lo = split(b, s).second;

   double e1 = res - a_hi * b_hi, 
      e2 = e1 - a_lo * b_hi, 
      e3 = e2 - a_hi * b_lo;

   roundoff = a_lo * b_lo - e3;
   return std::make_pair(res, roundoff);
}

void grow_expansion(std::vector<double> &e, double &b)
{
   double q = b;
   int n = e.size();
   for(int i = 0; i < n; i++)
   {
      std::pair<double, double> s = sum(q, e[i]);
      e[i] = s.second;
      q = s.first;
   }
   e.push_back(q);
}

std::vector<double> expand_sum(const std::vector<double> &e, const std::vector<double> &f)
{
   std::vector<double> res = e;
   int n = f.size();
   for(int i = 0; i < n; i++)
   {
      int m = res.size();
      double q = f[i];
      for(int j = i; j < m; j++)
      {
         std::pair<double, double> s = sum(q, res[j]);
         res[j] = s.second, q = s.first;
      }
      res.push_back(q);
   }
   return res;
}

//left turn & bounding box

bool bounding_box(double a, double b, double c, double d)
{
   if(a > b)
      std::swap(a, b);

   if (c > d)
      std::swap(c, d);

   if(a > d || c > b)
      return false;
   return true;

}

int left_turn(point a, point b, point c)
{
   double x = (b.x - a.x) * (c.y - a.y), y = (b.y - a.y) * (c.x - a.x);
   double eps = 4 * eps_m * (std::abs(x) + std::abs(y));
   double s = x - y;

   if(s > eps)
      return 1;
   if(s < -eps)
      return -1;

   interval_d inter_res = (interval_d(b.x) - interval_d(a.x)) * (interval_d(c.y) - interval_d(a.y)) -
      (interval_d(b.y) - interval_d(a.y)) * (interval_d(c.x) - interval_d(a.x));

   if(inter_res.lower() > 0)
      return 1;
   if(inter_res.upper() < 0)
      return -1;

   std::vector <std::pair<double, double>> mult;
   std::vector<double> e1, e2, e3;

   mult.push_back(multiply(b.x, c.y));
   mult.push_back(multiply(b.x, -a.y));
   mult.push_back(multiply(-a.x, c.y));
   mult.push_back(multiply(-b.y, c.x));
   mult.push_back(multiply(b.y, a.x));
   mult.push_back(multiply(a.y, c.x));

   std::vector<double> exp;
   for(int i = 0; i < 6; i++)
   {
      grow_expansion(exp, mult[i].first);
      grow_expansion(exp, mult[i].second);
   }

   for(int i = 11; i >= 0; i--)
   {
      if(exp[i] > 0)
         return 1;
      if(exp[i] < 0)
         return -1;
   }
   return 0;   

}

//intersection function
bool is_intersection(point a, point b, point c, point d)
{
   return bounding_box(a.x, b.x, c.x, d.x) && bounding_box(a.y, b.y, c.y, d.y) && left_turn(a, b, c) * left_turn(a, b, d) <= 0 && left_turn(c, d, a) * left_turn(c, d, b) <= 0;
}

#endif
